
# supress warnings (be careful while using this)
import warnings
warnings.filterwarnings('ignore')

# import common GIS tools
import numpy                   as np
import xarray                  as xr
import matplotlib.pyplot       as plt
import rasterio.features
import rioxarray               as rio
from matplotlib.cm import RdYlGn, Reds

# import Planetary Computer tools
import pystac_client
import planetary_computer as pc
import odc
from odc.stac import stac_load

# additional libraries
from datetime import date # date-related calculations

# GeoTiff images
import rasterio
from   osgeo import gdal

# data visualisation
from   matplotlib        import pyplot as plt
from   matplotlib.pyplot import figure
import matplotlib.image  as img
from   PIL               import Image

# model building
import ultralytics
from   ultralytics import YOLO
import labelme2yolo


# others
import os
import shutil
import zipfile


## Hurricane Maria - San Juan, Puerto Rico ##

## Defining the area of interest using latitude and longitude coordinates 
## at the center of the area we would like to observe

# Define the bounding box for the selected region 
min_lon = -66.19385887
min_lat =  18.27306794
max_lon = -66.069299
max_lat =  18.400288

# setting geographic boundary
bounds = (min_lon, min_lat, max_lon, max_lat)

# setting time window surrounding the storm landfall on September 20, 2017
time_window = "2017-04-08/2017-10-31"

# calculating days in time window
print(date(2017, 10, 31) - date(2017, 4, 8 ))

206 days, 0:00:00


## Using Planetary Computer's STAC catalog for items matching our query parameters above

# connecting to the planetary computer
stac = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")


# seaching for data
search = stac.search(collections = ["sentinel-2-l2a"],
                     bbox        = bounds,
                     datetime    = time_window)


# instantiating results list
items = list(search.get_all_items())


# summarizing results
print('This is the number of scenes that touch our region:',len(items))

This is the number of scenes that touch our region: 63


## Setting the resolution to 10 meters per pixel, and then convert this to degrees per pixel for the
## latitude-longitude coordinate system

# pixel resolution for the final product
resolution = 10  # meters per pixel 


# scaling to degrees per pizel
scale = resolution / 111320.0 # degrees per pixel for CRS:4326


## Use of Open Data Cube (ODC) for managing and analyzing geospatial data
## loading specific bands of data (like red, green, blue, near-infrared, and SCL)
## Filtering out clouds using the qa_pixel band and mapping at a resolution of 30 meters per pixel

xx = stac_load(
    items,
    bands      = ["red", "green", "blue", "nir", "SCL"],
    crs        = "EPSG:4326",                            # latitude-longitude
    resolution = scale,                                  # degrees
    chunks     = {"x": 2048, "y": 2048},
    dtype      = "uint16",
    patch_url  = pc.sign,
    bbox       = bounds
)


## Putting into comments since it will take some time to run

# # subsetting results for RGB
# plot_xx = xx[ ["red", "green", "blue"] ].to_array()


# # showing results 
# plot_xx.plot.imshow(col      = 'time', # time
#                     col_wrap = 4     , # four columns per row
#                     robust   = True  , 
#                     vmin     = 0     ,
#                     vmax     = 3000  )


# # rendering results
# plt.show()


# selecting a selected time slice to view a single RGB image and the cloud mask
time_slice = 29 # October 20, 2017 (post-storm)


# ## plotting an RGB real color image for a single date ##

# # setting plot size
# fig, ax = plt.subplots( figsize = (6, 10) )


# # preparing the plot
# xx.isel(time = time_slice)[ ["red", "green", "blue"] ].\
#     to_array().plot.imshow(robust = True,
#                            ax     = ax  ,
#                            vmin   = 0   ,
#                            vmax   = 3000)


# # titles and axis lables
# ax.set_title(label = f"RGB Color Results for Scene {time_slice}")
# ax.axis('off')


# # rendering results
# plt.show()

## This area of interest includes tropical forests and vegetation, some of which are located near Cupey and Caimito, 
## in the lower part of San Juan. These were significantly impacted as evidenced by the browning of vegetation. 
## According to NASA, nearly 60% of canopy trees in the region lost branches, snapped in half, or uprooted, with 
## trees that once had wide, spreading crowns now reduced to slender main trunks. Forests in Puerto Rico are now 
## approximately one-third shorter on average following Hurricane Maria (NASA, 2019).


# instantiating a colormap for SCL pixel classifications

scl_colormap = np.array(
    [
        [252,  40, 228, 255],  # 0  - NODATA - MAGENTA
        [255,   0,   4, 255],  # 1  - Saturated or Defective - RED
        [0  ,   0,   0, 255],  # 2  - Dark Areas - BLACK
        [97 ,  97,  97, 255],  # 3  - Cloud Shadow - DARK GREY
        [3  , 139,  80, 255],  # 4  - Vegetation - GREEN
        [192, 132,  12, 255],  # 5  - Bare Ground - BROWN
        [21 , 103, 141, 255],  # 6  - Water - BLUE
        [117,   0,  27, 255],  # 7  - Unclassified - MAROON
        [208, 208, 208, 255],  # 8  - Cloud - LIGHT GREY
        [244, 244, 244, 255],  # 9  - Definitely Cloud - WHITE
        [195, 231, 240, 255],  # 10 - Thin Cloud - LIGHT BLUE
        [222, 157, 204, 255],  # 11 - Snow or Ice - PINK
    ],
    dtype="uint8",
)


# function for color encoding
def colorize(xx, colormap):
    return xr.DataArray( colormap[xx.data],
                         coords = xx.coords,
                         dims   = (*xx.dims, "band") )


# filtering out water, etc.
filter_values = [0, 1, 3, 6, 8, 9, 10]

cloud_mask = ~xx.SCL.isin(filter_values) # this means not in filter_values


# appling cloud mask (filtering out clouds, cloud shadows, and water)

# storing as 16-bit integers
cleaned_data = xx.where(cloud_mask).astype("uint16")


# # converting SCL to RGB
# scl_rgba_clean = colorize(xx       = cleaned_data.isel(time = time_slice).SCL.compute(), 
#                           colormap = scl_colormap)


# # setting figure size
# plt.figure(figsize = (6, 10))
# plt.imshow(scl_rgba_clean)


# # titles and axis labels
# plt.title(label = "Cloud / Shadows / Water Mask (MAGENTA)")
# plt.axis('off')


# # rendering the plot
# plt.show()

## Code above creates a colormap for different land cover types based on Sentinel-2 Scene Classification Layer (SCL) 
## values. It then applies the colormap to the data, highlighting vegetation, bare ground, and water in distinct colors. 
## It filters out clouds, cloud shadows, and no data from the image represented by the magenta color, displaying the remaining land cover types in the area of interest.


# preparing two time steps compare NDVI outputs
first_time  = 0  # April 08, 2017
second_time = 29 # October 20, 2017


# # Plots of NDVI at two different time slices
# # Clouds, Shadows and Water are masked

# # setting figure size
# fig, ax = plt.subplots(1, 2, figsize = (15, 10))


# # first image data
# ndvi_image = (cleaned_data.nir - cleaned_data.red) / (cleaned_data.nir + cleaned_data.red)
# ndvi_image.isel(time = first_time ).plot(ax = ax[0],
#                                          vmin = 0.0,
#                                          vmax = 0.8,
#                                          cmap = "RdYlGn")


# # second image data
# ndvi_image.isel(time = second_time).plot(ax = ax[1],
#                                          vmin = 0.0,
#                                          vmax = 0.8,
#                                          cmap = "RdYlGn")


# # axis labels
# ax[0].set_title(label = 'NDVI-Time #1'), ax[0].xaxis.set_visible(False), ax[0].yaxis.set_visible(False)
# ax[1].set_title(label = 'NDVI-Time #2'), ax[1].xaxis.set_visible(False), ax[1].yaxis.set_visible(False)


# # rendering the plot
# plt.show()

## The figure is from April 08 and October 20, 2017. It creates plots of NDVI values for selected time window,
## masking out clouds, shadows, and no data. Values indicate the health and density of vegetation, with higher 
## values typically indicating healthier vegetation (dark green). Moreover, it illustrates how vegetation cover 
## changed post-storm, evident in the shift to lighter green, red, and orange shades, providing insights into the 
## sudden drop in NDVI after the storm.


# function for calculating NDVI anomalies
def NDVI(dataset):
    return (dataset.nir - dataset.red) / (dataset.nir + dataset.red)


# running comparison
ndvi_clean = NDVI(cleaned_data)

# calculating difference
ndvi_pre     = ndvi_clean.isel(time = first_time)
ndvi_post    = ndvi_clean.isel(time = second_time)
ndvi_anomaly = ndvi_post - ndvi_pre


# all areas of water or clouds will be black
RdYlGn.set_bad('black',1.)


# reversing the colormap for reds
Reds_reverse = "Reds_r"


## plotting NDVI anomaly
plt.figure( figsize = (6,10) )
ndvi_anomaly.plot(vmin = -0.2, vmax=0.0, cmap = Reds_reverse, add_colorbar=False)


# titles and labels
plt.title (label  = "NDVI Anomaly")
plt.xlabel(xlabel = "Longitude")
plt.ylabel(ylabel = "Latitude")

# Save the current figure as a variable
ndvi_plot_01 = plt.gcf()

# Display the plot
plt.show()


# Define the bounding box for the selected region
min_lon = -66.19385887
min_lat =  18.27306794
max_lon = -66.08007533
max_lat =  18.48024350

# setting geographic boundary
bounds = (min_lon, min_lat, max_lon, max_lat)


# setting time window surrounding the storm landfall on September 20, 2017
time_window = "2017-04-08/2017-11-30"

# calculating days in time window
print(date(2017, 11, 1) - date(2017, 4, 1 ))

214 days, 0:00:00


## Using Planetary Computer's STAC catalog for items matching our query parameters above

# connecting to the planetary computer
stac = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")


# seaching for data
search = stac.search(collections = ["sentinel-2-l2a"],
                     bbox        = bounds,
                     datetime    = time_window)


# instantiating results list
items = list(search.get_all_items())


# summarizing results
print('This is the number of scenes that touch our region:',len(items))

This is the number of scenes that touch our region: 73


## Setting the resolution to 10 meters per pixel, and then convert this to degrees per pixel for the
## latitude-longitude coordinate system

# pixel resolution for the final product
resolution = 10  # meters per pixel 


# scaling to degrees per pizel
scale = resolution / 111320.0 # degrees per pixel for CRS:4326


## Use of Open Data Cube (ODC) for managing and analyzing geospatial data
## loading specific bands of data (like red, green, blue, near-infrared, and SCL)
## Filtering out clouds using the qa_pixel band and mapping at a resolution of 30 meters per pixel

xx = stac_load(
    items,
    bands      = ["red", "green", "blue", "nir", "SCL"],
    crs        = "EPSG:4326",                            # latitude-longitude
    resolution = scale,                                  # degrees
    chunks     = {"x": 2048, "y": 2048},
    dtype      = "uint16",
    patch_url  = pc.sign,
    bbox       = bounds
)


## Putting into comments since it will take some time to run

# # subsetting results for RGB
# plot_xx = xx[ ["red", "green", "blue"] ].to_array()


# # showing results 
# plot_xx.plot.imshow(col      = 'time', # time
#                     col_wrap = 4     , # four columns per row
#                     robust   = True  , 
#                     vmin     = 0     ,
#                     vmax     = 3000  )


# # rendering results
# plt.show()

## This visual focuses in the middle part where houses are located near forested areas. Geographic area and 
## time frame, relevant satellite data, and various angles were displayed. Visualizations have more clouds
## and not a good angle, usually either cropped or not captured at all. Zooming in on the area was done to 
## assess the extent of damage. Same pre- and post-dates were chosen for consistency in assessing damage only 
## varying the coordinates.


# selecting a selected time slice to view a single RGB image and the cloud mask
time_slice = 30 # October 20, 2017 (post storm)


# ## plotting an RGB real color image for a single date ##

# # setting plot size
# fig, ax = plt.subplots( figsize = (6, 10) )


# # preparing the plot
# xx.isel(time = time_slice)[ ["red", "green", "blue"] ].\
#     to_array().plot.imshow(robust = True,
#                            ax     = ax  ,
#                            vmin   = 0   ,
#                            vmax   = 3000)


# # titles and axis lables
# ax.set_title(label = f"RGB Color Results for Scene {time_slice}")
# ax.axis('off')


# # rendering results
# plt.show()

## The highlighted area of interest includes tropical forests, vegetation, and surrounding suburbs, indicated by 
## white tones. Notably, there is observable browning of vegetation near forested regions, suburbs, and infrastructure 
## along the coastal areas. These findings align with Hosannah, N. et al. (2019), who noted a significant increase in 
## barren and pasture land cover post-storm, with forest cover being converted to bare soil (13%–41%) or pasture in 
## both forested areas and suburbs.


# filtering out water, etc.
filter_values = [0, 1, 3, 6, 8, 9, 10]

cloud_mask = ~xx.SCL.isin(filter_values) # this means not in filter_values


# appling cloud mask (filtering out clouds, cloud shadows, and water)

# storing as 16-bit integers
cleaned_data = xx.where(cloud_mask).astype("uint16")


## Putting into comments; same code as Viz 1

# # converting SCL to RGB
# scl_rgba_clean = colorize(xx       = cleaned_data.isel(time = time_slice).SCL.compute(), 
#                           colormap = scl_colormap)


# # setting figure size
# plt.figure(figsize = (6, 10))
# plt.imshow(scl_rgba_clean)


# # titles and axis labels
# plt.title(label = "Cloud / Shadows / Water Mask (MAGENTA)")
# plt.axis('off')


# # rendering the plot
# plt.show()


# preparing two time steps compare NDVI outputs
first_time  = 0  # April 08, 2017
second_time = 30 # October 20, 2017


## Putting into comments; same code as Viz 1

# # Plots of NDVI at two different time slices
# # Clouds, Shadows and Water are masked

# # setting figure size
# fig, ax = plt.subplots(1, 2, figsize = (15, 10))


# # first image data
# ndvi_image = (cleaned_data.nir - cleaned_data.red) / (cleaned_data.nir + cleaned_data.red)
# ndvi_image.isel(time = first_time ).plot(ax = ax[0],
#                                          vmin = 0.0,
#                                          vmax = 0.8,
#                                          cmap = "RdYlGn")


# # second image data
# ndvi_image.isel(time = second_time).plot(ax = ax[1],
#                                          vmin = 0.0,
#                                          vmax = 0.8,
#                                          cmap = "RdYlGn")


# # axis labels
# ax[0].set_title(label = 'NDVI-Time #1'), ax[0].xaxis.set_visible(False), ax[0].yaxis.set_visible(False)
# ax[1].set_title(label = 'NDVI-Time #2'), ax[1].xaxis.set_visible(False), ax[1].yaxis.set_visible(False)


# # rendering the plot
# plt.show()


# running comparison
ndvi_clean = NDVI(cleaned_data)

# calculating difference
ndvi_pre     = ndvi_clean.isel(time = first_time)
ndvi_post    = ndvi_clean.isel(time = second_time)
ndvi_anomaly = ndvi_post - ndvi_pre


# all areas of water or clouds will be black
RdYlGn.set_bad('black',1.)


# reversing the colormap for reds
Reds_reverse = "Reds_r"


## plotting NDVI anomaly
plt.figure( figsize = (6,10) )
ndvi_anomaly.plot(vmin = -0.2, vmax=0.0, cmap = Reds_reverse, add_colorbar=False)


# titles and labels
plt.title (label  = "NDVI Anomaly")
plt.xlabel(xlabel = "Longitude")
plt.ylabel(ylabel = "Latitude")

# Save the current figure as a variable
ndvi_plot_02 = plt.gcf()

# Display the plot
plt.show()


# Define the bounding box for the selected region
min_lon = -66.1327369
min_lat =  18.4063942
max_lon = -66.0673603
max_lat =  18.4784524


# setting geographic boundary
bounds = (min_lon, min_lat, max_lon, max_lat)


# setting time window surrounding the storm landfall on September 20, 2017
time_window = "2017-04-01/2017-11-01"

# calculating days in time window
print(date(2017, 11, 1) - date(2017, 4, 1 ))

214 days, 0:00:00


## Using Planetary Computer's STAC catalog for items matching our query parameters above

# connecting to the planetary computer
stac = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")


# seaching for data
search = stac.search(collections = ["sentinel-2-l2a"],
                     bbox        = bounds,
                     datetime    = time_window)


# instantiating results list
items = list(search.get_all_items())


# summarizing results
print('This is the number of scenes that touch our region:',len(items))

This is the number of scenes that touch our region: 45


## Setting the resolution to 10 meters per pixel, and then convert this to degrees per pixel for the
## latitude-longitude coordinate system

# pixel resolution for the final product
resolution = 10  # meters per pixel 


# scaling to degrees per pizel
scale = resolution / 111320.0 # degrees per pixel for CRS:4326


## Use of Open Data Cube (ODC) for managing and analyzing geospatial data
## loading specific bands of data (like red, green, blue, near-infrared, and SCL)
## Filtering out clouds using the qa_pixel band and mapping at a resolution of 30 meters per pixel

xx = stac_load(
    items,
    bands      = ["red", "green", "blue", "nir", "SCL"],
    crs        = "EPSG:4326",                            # latitude-longitude
    resolution = scale,                                  # degrees
    chunks     = {"x": 2048, "y": 2048},
    dtype      = "uint16",
    patch_url  = pc.sign,
    bbox       = bounds
)


## Putting into comments since it will take some time to run

# # subsetting results for RGB
# plot_xx = xx[ ["red", "green", "blue"] ].to_array()


# # showing results 
# plot_xx.plot.imshow(col      = 'time', # time
#                     col_wrap = 4     , # four columns per row
#                     robust   = True  , 
#                     vmin     = 0     ,
#                     vmax     = 3000  )


# # rendering results
# plt.show()

## This visualization targets an area in northern San Juan, Puerto Rico, featuring critical infrastructure near 
## coastal regions. It defines a specific geographical area and time range and collects pertinent satellite data.
## However, many visualizations are obscured by clouds, and some lack data entirely. To evaluate damage, we zoomed 
## into the area, maintaining consistent pre- and post-dates while varying coordinates to assess damage extent in
## this area.


# selecting a selected time slice to view a single RGB image and the cloud mask
time_slice = 21 # October 20, 2017 (post storm)


# ## plotting an RGB real color image for a single date ##

# # setting plot size
# fig, ax = plt.subplots( figsize = (6, 10) )


# # preparing the plot
# xx.isel(time = time_slice)[ ["red", "green", "blue"] ].\
#     to_array().plot.imshow(robust = True,
                           
#                            ax     = ax  ,
#                            vmin   = 0   ,
#                            vmax   = 3000)


# # titles and axis lables
# ax.set_title(label = f"RGB Color Results for Scene {time_slice}")
# ax.axis('off')


# # rendering results
# plt.show()

## It is evident that vegetation near coastal infrastructures has also shown signs of browning. This indicates 
## potential damage to houses and buildings. This image can help identify the specific types of buildings, 
## critical infrastructures, and areas most affected by the storm. Moreover, beach erosion and its significant 
## impact on coastal vulnerability can be assessed, as studied by Valentin & Muller (2020), who found an average 
## beach loss of 3 to 5 meters along the hurricane's path and up to 40 meters at certain beaches


# filtering out water, etc.
filter_values = [0, 1, 3, 6, 8, 9, 10]

cloud_mask = ~xx.SCL.isin(filter_values) # this means not in filter_values


# appling cloud mask (filtering out clouds, cloud shadows, and water)

# storing as 16-bit integers
cleaned_data = xx.where(cloud_mask).astype("uint16")


## Putting into comments; same code as Viz 1/2

# # converting SCL to RGB
# scl_rgba_clean = colorize(xx       = cleaned_data.isel(time = time_slice).SCL.compute(), 
#                           colormap = scl_colormap)


# # setting figure size
# plt.figure(figsize = (6, 10))
# plt.imshow(scl_rgba_clean)


# # titles and axis labels
# plt.title(label = "Cloud / Shadows / Water Mask (MAGENTA)")
# plt.axis('off')


# # rendering the plot
# plt.show()


# preparing two time steps compare NDVI outputs
first_time  = 0  # April 08, 2017
second_time = 21 # October 20, 2017


## Putting into comments; same code as Viz 1/2

# # Plots of NDVI at two different time slices
# # Clouds, Shadows and Water are masked

# # setting figure size
# fig, ax = plt.subplots(1, 2, figsize = (15, 10))


# # first image data
# ndvi_image = (cleaned_data.nir - cleaned_data.red) / (cleaned_data.nir + cleaned_data.red)
# ndvi_image.isel(time = first_time ).plot(ax = ax[0],
#                                          vmin = 0.0,
#                                          vmax = 0.8,
#                                          cmap = "RdYlGn")


# # second image data
# ndvi_image.isel(time = second_time).plot(ax = ax[1],
#                                          vmin = 0.0,
#                                          vmax = 0.8,
#                                          cmap = "RdYlGn")


# # axis labels
# ax[0].set_title(label = 'NDVI-Time #1'), ax[0].xaxis.set_visible(False), ax[0].yaxis.set_visible(False)
# ax[1].set_title(label = 'NDVI-Time #2'), ax[1].xaxis.set_visible(False), ax[1].yaxis.set_visible(False)


# # rendering the plot
# plt.show()

## Using the same time window as before, this visualization generates plots of NDVI values while excluding clouds, 
## shadows, water, and dark areas near the coastal area. NDVI values indicate vegetation health and density, with 
## higher values suggesting healthier vegetation (dark green). The plots demonstrate changes in vegetation cover 
## following the storm near Guaynabo, Puerto Nuevo, and Catano, among others, highlighting shifts to lighter green, 
## red, yellow, and orange shades in forested, suburban, and coastal regions. These changes provide valuable insights 
## into the sudden decline in NDVI post-storm.


# running comparison
ndvi_clean = NDVI(cleaned_data)

# calculating difference
ndvi_pre     = ndvi_clean.isel(time = first_time)
ndvi_post    = ndvi_clean.isel(time = second_time)
ndvi_anomaly = ndvi_post - ndvi_pre


# all areas of water or clouds will be black
RdYlGn.set_bad('black',1.)


# reversing the colormap for reds
Reds_reverse = "Reds_r"


## plotting NDVI anomaly
plt.figure( figsize = (6,10) )
ndvi_anomaly.plot(vmin = -0.2, vmax=0.0, cmap = Reds_reverse, add_colorbar=False)


# titles and labels
plt.title (label  = "NDVI Anomaly")
plt.xlabel(xlabel = "Longitude")
plt.ylabel(ylabel = "Latitude")

# Save the current figure as a variable
ndvi_plot_03 = plt.gcf()

# Display the plot
plt.show()


#!labelme2yolo --json_dir /path/to/labelme_json_dir/
!labelme2yolo --json_dir "./labelme_json_dir/"

Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--
Converting... -- -------------------------------------   7% 0:00:01
Converting... ----- ----------------------------------  15% 0:00:02
Converting... -------- -------------------------------  21% 0:00:02
Converting... ---------- -----------------------------  26% 0:00:02
Converting... -------------- -------------------------  35% 0:00:01
Converting... ----------------- ----------------------  45% 0:00:01
Converting... -------------------- -------------------  52% 0:00:01
Converting... ----------------------- ----------------  59% 0:00:01
Converting... --------------------------- ------------  70% 0:00:01
Converting... ------------------------------- --------  79% 0:00:01
Converting... ---------------------------------- -----  86% 0:00:01
Converting... ------------------------------------- --  93% 0:00:01
Converting... ---------------------------------------- 100% 0:00:00

Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--
Converting... ----------- ----------------------------  29% 0:00:01
Converting... --------------------- ------------------  54% 0:00:01
Converting... ------------------------------------ ---  92% 0:00:01
Converting... ---------------------------------------- 100% 0:00:00

Converting... ----------------------------------------   0% -:--:--
Converting... ----------------------------------------   0% -:--:--

INFO:labelme2yolo:Searching label list from json files ...
INFO:labelme2yolo:Converting train set ...
INFO:labelme2yolo:Converting val set ...
INFO:labelme2yolo:Converting test set ...


#Loading the model
model = YOLO('yolov8n.pt')

# Display model information (optional)
model.info()

YOLOv8n summary: 225 layers, 3157200 parameters, 0 gradients, 8.9 GFLOPs

(225, 3157200, 0, 8.8575488)


# Train the model on the dataset for 50 epochs
results = model.train(data   = "./datasets/dataset.yaml", 
                      epochs = 50  ,
                      imgsz  = 512)

New https://pypi.org/project/ultralytics/8.2.4 available  Update with 'pip install -U ultralytics'
Ultralytics YOLOv8.1.47  Python-3.11.5 torch-2.2.2+cpu CPU (13th Gen Intel Core(TM) i7-1355U)
engine\trainer: task=detect, mode=train, model=yolov8n.pt, data=./datasets/dataset.yaml, epochs=50, time=None, patience=100, batch=16, imgsz=512, save=True, save_period=-1, cache=False, device=None, workers=0, project=None, name=train642, exist_ok=False, pretrained=True, optimizer=auto, verbose=True, seed=0, deterministic=True, single_cls=False, rect=False, cos_lr=False, close_mosaic=10, resume=False, amp=True, fraction=1.0, profile=False, freeze=None, multi_scale=False, overlap_mask=True, mask_ratio=4, dropout=0.0, val=True, split=val, save_json=False, save_hybrid=False, conf=None, iou=0.7, max_det=300, half=False, dnn=False, plots=True, source=None, vid_stride=1, stream_buffer=False, visualize=False, augment=False, agnostic_nms=False, classes=None, retina_masks=False, embed=None, show=False, save_frames=False, save_txt=False, save_conf=False, save_crop=False, show_labels=True, show_conf=True, show_boxes=True, line_width=None, format=torchscript, keras=False, optimize=False, int8=False, dynamic=False, simplify=False, opset=None, workspace=4, nms=False, lr0=0.01, lrf=0.01, momentum=0.937, weight_decay=0.0005, warmup_epochs=3.0, warmup_momentum=0.8, warmup_bias_lr=0.0, box=7.5, cls=0.5, dfl=1.5, pose=12.0, kobj=1.0, label_smoothing=0.0, nbs=64, hsv_h=0.015, hsv_s=0.7, hsv_v=0.4, degrees=0.0, translate=0.1, scale=0.5, shear=0.0, perspective=0.0, flipud=0.0, fliplr=0.5, bgr=0.0, mosaic=0.0, mixup=0.0, copy_paste=0.0, auto_augment=randaugment, erasing=0.4, crop_fraction=1.0, cfg=None, tracker=botsort.yaml, save_dir=runs\detect\train642

                   from  n    params  module                                       arguments                     
  0                  -1  1       464  ultralytics.nn.modules.conv.Conv             [3, 16, 3, 2]                 
  1                  -1  1      4672  ultralytics.nn.modules.conv.Conv             [16, 32, 3, 2]                
  2                  -1  1      7360  ultralytics.nn.modules.block.C2f             [32, 32, 1, True]             
  3                  -1  1     18560  ultralytics.nn.modules.conv.Conv             [32, 64, 3, 2]                
  4                  -1  2     49664  ultralytics.nn.modules.block.C2f             [64, 64, 2, True]             
  5                  -1  1     73984  ultralytics.nn.modules.conv.Conv             [64, 128, 3, 2]               
  6                  -1  2    197632  ultralytics.nn.modules.block.C2f             [128, 128, 2, True]           
  7                  -1  1    295424  ultralytics.nn.modules.conv.Conv             [128, 256, 3, 2]              
  8                  -1  1    460288  ultralytics.nn.modules.block.C2f             [256, 256, 1, True]           
  9                  -1  1    164608  ultralytics.nn.modules.block.SPPF            [256, 256, 5]                 
 10                  -1  1         0  torch.nn.modules.upsampling.Upsample         [None, 2, 'nearest']          
 11             [-1, 6]  1         0  ultralytics.nn.modules.conv.Concat           [1]                           
 12                  -1  1    148224  ultralytics.nn.modules.block.C2f             [384, 128, 1]                 
 13                  -1  1         0  torch.nn.modules.upsampling.Upsample         [None, 2, 'nearest']          
 14             [-1, 4]  1         0  ultralytics.nn.modules.conv.Concat           [1]                           
 15                  -1  1     37248  ultralytics.nn.modules.block.C2f             [192, 64, 1]                  
 16                  -1  1     36992  ultralytics.nn.modules.conv.Conv             [64, 64, 3, 2]                
 17            [-1, 12]  1         0  ultralytics.nn.modules.conv.Concat           [1]                           
 18                  -1  1    123648  ultralytics.nn.modules.block.C2f             [192, 128, 1]                 
 19                  -1  1    147712  ultralytics.nn.modules.conv.Conv             [128, 128, 3, 2]              
 20             [-1, 9]  1         0  ultralytics.nn.modules.conv.Concat           [1]                           
 21                  -1  1    493056  ultralytics.nn.modules.block.C2f             [384, 256, 1]                 
 22        [15, 18, 21]  1    752092  ultralytics.nn.modules.head.Detect           [4, [64, 128, 256]]           
Model summary: 225 layers, 3011628 parameters, 3011612 gradients, 8.2 GFLOPs

Transferred 355/355 items from pretrained weights
Freezing layer 'model.22.dfl.conv.weight'

train: Scanning C:\Users\joann\Documents\00 Hult 2023 2024\Spring 1 2024 Courses\Business Challenge III\From Chase\data

train: New cache created: C:\Users\joann\Documents\00 Hult 2023 2024\Spring 1 2024 Courses\Business Challenge III\From Chase\datasets\train.cache
WARNING  Box and segment counts should be equal, but got len(segments) = 680, len(boxes) = 15337. To resolve this only boxes will be used and all segments will be removed. To avoid this please supply either a detect or segment dataset, not a detect-segment mixed dataset.

val: Scanning C:\Users\joann\Documents\00 Hult 2023 2024\Spring 1 2024 Courses\Business Challenge III\From Chase\datase

val: New cache created: C:\Users\joann\Documents\00 Hult 2023 2024\Spring 1 2024 Courses\Business Challenge III\From Chase\datasets\val.cache
WARNING  Box and segment counts should be equal, but got len(segments) = 163, len(boxes) = 3303. To resolve this only boxes will be used and all segments will be removed. To avoid this please supply either a detect or segment dataset, not a detect-segment mixed dataset.
Plotting labels to runs\detect\train642\labels.jpg... 
optimizer: 'optimizer=auto' found, ignoring 'lr0=0.01' and 'momentum=0.937' and determining best 'optimizer', 'lr0' and 'momentum' automatically... 
optimizer: AdamW(lr=0.00125, momentum=0.9) with parameter groups 57 weight(decay=0.0), 64 weight(decay=0.0005), 63 bias(decay=0.0)
Image sizes 512 train, 512 val
Using 0 dataloader workers
Logging results to runs\detect\train642
Starting training for 50 epochs...

      Epoch    GPU_mem   box_loss   cls_loss   dfl_loss  Instances       Size

       1/50         0G      0.767     0.8557     0.9686        644        512: 100%|██████████| 26/26 [03:06<00:00,  7.
                 Class     Images  Instances      Box(P          R      mAP50  mAP50-95): 100%|██████████| 4/4 [00:24<0

                   all        104       3303      0.669      0.698      0.725      0.563

      Epoch    GPU_mem   box_loss   cls_loss   dfl_loss  Instances       Size

       2/50         0G     0.8003     0.8974     0.9787        432        512: 100%|██████████| 26/26 [03:01<00:00,  6.
                 Class     Images  Instances      Box(P          R      mAP50  mAP50-95): 100%|██████████| 4/4 [00:23<0


%matplotlib inline

figure(figsize=(15, 10), dpi=80)
# reading the image 
results = img.imread('./runs/detect/train642/results.png')

# displaying the image 
plt.imshow(results) 
plt.show()


figure(figsize=(20,15), dpi=80)  

# reading the image 
cf = img.imread('./runs/detect/train642/confusion_matrix.png') 

# displaying the image 
plt.imshow(cf) 
plt.show()


# Confusion Matrix for Class 1 (2x2)
confusion_matrix_one = np.array([ [59, 104],
                                  [73,  4429] ] )

# Select rows and columns corresponding to class 1
class_index = 0
tn_1 = confusion_matrix_one[class_index, class_index]
fp_1 = confusion_matrix_one[:, class_index].sum() - tp_1
fn_1 = confusion_matrix_one[class_index, :].sum() - tp_1
tp_1 = confusion_matrix_one.sum() - tn_1 - fp_1 - fn_1

# Construct one-vs-all confusion matrix for class 1
class_one = np.array([[tp_1, fp_1],
                      [fn_1, tn_1]])

print("Confusion Matrix for Class 1 (undamaged residential building)")

# Unpacking the one-vs-all confusion matrix for Class 1
tn_1, fp_1, fn_1, tp_1 = class_one.ravel()

# Printing each result one-by-one
print(f"""
True Positives: {tp_1}
False Positives: {fp_1}
True Negatives: {tn_1}
False Negatives: {fn_1}
""")

# Calculate precision and recall for Class 1 (undamaged residential building)
precision_1 = tp_1 / (tp_1 + fp_1)
recall_1 = tp_1 / (tp_1 + fn_1)

print(f"Precision for Class 1: {precision_1:.4f}")
print(f"Recall for Class 1: {recall_1:.4f}")

Confusion Matrix for Class 1 (undamaged residential building)

True Positives: 59
False Positives: 73
True Negatives: 4429
False Negatives: 104

Precision for Class 1: 0.4470
Recall for Class 1: 0.3620


# Confusion Matrix for Class 2 (2x2)
confusion_matrix_two = np.array([ [379, 473],
                                  [173,  3640] ] )

# Select rows and columns corresponding to class 2
class_index = 0
tn_2 = confusion_matrix_two[class_index, class_index]
fp_2 = confusion_matrix_two[:, class_index].sum() - tp_2
fn_2 = confusion_matrix_two[class_index, :].sum() - tp_2
tp_2 = confusion_matrix_two.sum() - tn_2 - fp_2 - fn_2

# Construct one-vs-all confusion matrix for class 2
class_two = np.array([[tp_2, fp_2],
                      [fn_2, tn_2]])

print("Confusion Matrix for Class 2 (damaged residential building)")

# Unpacking the one-vs-all confusion matrix for Class 2
tn_2, fp_2, fn_2, tp_2 = class_two.ravel()

# Printing each result one-by-one
print(f"""
True Positives: {tp_2}
False Positives: {fp_2}
True Negatives: {tn_2}
False Negatives: {fn_2}
""")

# Calculate precision and recall for Class 2 (damaged residential buildings)
precision_2 = tp_2 / (tp_2 + fp_2)
recall_2 = tp_2 / (tp_2 + fn_2)

print(f"Precision for Class 2: {precision_2:.4f}")
print(f"Recall for Class 2: {recall_2:.4f}")

Confusion Matrix for Class 2 (damaged residential building)

True Positives: 379
False Positives: 173
True Negatives: 3640
False Negatives: 473

Precision for Class 2: 0.6866
Recall for Class 2: 0.4448


# Confusion Matrix for Class 3 (2x2)
confusion_matrix_three = np.array([ [173, 388],
                                    [134, 3970] ] )

# Select rows and columns corresponding to class 3
class_index = 0
tn_3 = confusion_matrix_three[class_index, class_index]
fp_3 = confusion_matrix_three[:, class_index].sum() - tp_3
fn_3 = confusion_matrix_three[class_index, :].sum() - tp_3
tp_3 = confusion_matrix_three.sum() - tn_3 - fp_3 - fn_3

# Construct one-vs-all confusion matrix for class 3
class_three = np.array([[tp_3, fp_3],
                      [fn_3, tn_3]])

print("Confusion Matrix for Class 3 (undamaged commercial building)")

# Unpacking the one-vs-all confusion matrix for Class 3
tn_3, fp_3, fn_3, tp_3 = class_three.ravel()

# Printing each result one-by-one
print(f"""
True Positives: {tp_3}
False Positives: {fp_3}
True Negatives: {tn_3}
False Negatives: {fn_3}
""")

# Calculate precision and recall for Class 3
precision_3 = tp_3 / (tp_3 + fp_3)
recall_3 = tp_3 / (tp_3 + fn_3)

print(f"Precision for Class 3: {precision_3:.4f}")
print(f"Recall for Class 3: {recall_3:.4f}")

Confusion Matrix for Class 3 (undamaged commercial building)

True Positives: 173
False Positives: 134
True Negatives: 3970
False Negatives: 388

Precision for Class 3: 0.5635
Recall for Class 3: 0.3084


# Confusion Matrix for Class 4 (2x2)
confusion_matrix_four = np.array([ [1853, 810],
                                    [459, 1543] ] )

# Select rows and columns corresponding to class 4
class_index = 0
tn_4 = confusion_matrix_four[class_index, class_index]
fp_4 = confusion_matrix_four[:, class_index].sum() - tp_4
fn_4 = confusion_matrix_four[class_index, :].sum() - tp_4
tp_4 = confusion_matrix_four.sum() - tn_4 - fp_4 - fn_4

# Construct one-vs-all confusion matrix for class 4
class_four = np.array([[tp_4, fp_4],
                      [fn_4, tn_4]])

print("Confusion Matrix for Class 4 (damaged commercial building)")

# Unpacking the one-vs-all confusion matrix for Class 4
tn_4, fp_4, fn_4, tp_4 = class_four.ravel()

# Printing each result one-by-one
print(f"""
True Positives: {tp_4}
False Positives: {fp_4}
True Negatives: {tn_4}
False Negatives: {fn_4}
""")

# Calculate precision and recall for Class 4
precision_4 = tp_4 / (tp_4 + fp_4)
recall_4 = tp_4 / (tp_4 + fn_4)

print(f"Precision for Class 2: {precision_4:.4f}")
print(f"Recall for Class 2: {recall_4:.4f}")

Confusion Matrix for Class 4 (damaged commercial building)

True Positives: 1853
False Positives: 459
True Negatives: 1543
False Negatives: 810

Precision for Class 2: 0.8015
Recall for Class 2: 0.6958

Tropical Cyclone Damage Assessment

BCH-7815 | Business Challenge III

Introduction

Top Three Actionable Insights

Libraries and packages

Part I. Pre- and Post-Event NDVI Analysis

Visualization 01

Accessing Satellite Data

Viewing RGB (real color) images from the time series

Applying Cloud Filtering and Masking

Normalized Difference Vegetation Index (NDVI)

NDVI Change Product

Analysis (Visualization 01)

Visualization 02

Accessing Satellite Data

Viewing RGB (real color) images from the time series

Applying Cloud Filtering and Masking

Normalized Difference Vegetation Index (NDVI)

NDVI Change Product

Analysis (Visualization 02)

Visualization 03

Accessing Satellite Data

Viewing RGB (real color) images from the time series

Applying Cloud Filtering and Masking

Normalized Difference Vegetation Index (NDVI)

NDVI Change Product

Analysis (Visualization 03)

Part II. Object Detection and Model Building

Generating Training and Testing Data along with the Configuration File

Model Building and Model Training

Model Evaluation

Part III Model Analysis

Confusion Matrix

Converting Confusion Matrix to One Class

Class 1 Undamaged residential buildings

Class 2 Damaged residential buildings

Class 3 Undamaged commercial building

Class 4 Damaged commercial building

Confusion Matrix Analysis

Conclusion

Steps your team would implement/improve if you were given three months to work on this project

Part IV. Feedback to EY

References